      subroutine parmovgc 
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use corgan_com_M 
      use cindex_com_M 
      use numpar_com_M 
      use cprplt_com_M, ONLY:  
      use cophys_com_M 
      use ctemp_com_M, ONLY: fmf, fef, fvf, fbxf, fbyf, fbzf, twx 
      use blcom_com_M, ONLY: drift, qom, t_wall, iphd2, wate, divpix, bxv, &
         byv, bzv, gradxb, gradyb, gradzb, ex, ey, ez, x, y, z, area_x, &
         area_y, area_z, tsix, tsiy, tsiz, iphead, ijkcell, ijktmp2, ijktmp3, &
         ijktmp4, ijkctmp, nux, nuy, nuz, etax, etay, etaz, px, py, pz, pmu, &
         pxi, peta, pzta, up, vp, wp, link, ico, xptilde, yptilde, zptilde, &
         uptilde, vptilde, wptilde, bxpn, bypn, bzpn, expn, eypn, ezpn, vrms 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
!
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer , dimension(8) :: iv 
      integer , dimension(4) :: lioinput, lioprint 
      integer :: lpx, n, np, is, inew, jnew, knew, ijknew 
      real(double), dimension(1) :: ixi1, ixi2, ieta1, ieta2 
      real(double) :: ixi4, ieta4, ixi5, ieta5 
      real(double), dimension(8) :: wght 
      real(double), dimension(100) :: fpxf, fpyf, fpzf, fpxft, fpyft, fpzft, &
         ucm, vcm, wcm, cm 
      real(double), dimension(3,20) :: wsin, wcos 
      real(double), dimension(3,12,20) :: tsi, rot 
      real(double) :: nuxp, nuyp, nuzp, nu, dto2, qomdt, rqom, omdtsq, denom, &
         ehatx, ehaty, ehatz, rbsq, udotb 
      logical :: expnd, nomore 
      character :: name*80 
!-----------------------------------------------
!
!
!
!     *******************************************
!
!     solves guiding-center equations
!     a modified leap-frog algorithm is used
!
!     *********************************************
      lpx = nphist + 1 
      dto2 = 0.5*dt 
!
!dir$ ivdep
!
      iphd2(ijkcell(:ncells)) = 0 
!
!
    1 continue 
      nomore = .TRUE. 
!
!     check for more particles to process
!
      ncellsp = 0 
      do n = 1, ncells 
  211    continue 
         if (iphead(ijkcell(n)) <= 0) cycle  
         ijk = ijkcell(n) 
         np = iphead(ijk) 
         is = ico(np) 
         if (drift(is)) then 
!
            ncellsp = ncellsp + 1 
            ijkctmp(ncellsp) = ijkcell(n) 
!
         else 
!
            iphead(ijk) = link(np) 
            link(np) = iphd2(ijk) 
            iphd2(ijk) = np 
!
            go to 211 
!
         endif 
      end do 
!
!     if there are more, process particles as though there were one
!     in every cell in the mesh
!
      if (ncellsp /= 0) then 
!
         nomore = .FALSE. 
!
!
!
         do n = 1, ncellsp 
!
            ijk = ijkctmp(n) 
            np = iphead(ijk) 
!
            px(np) = px(np) + up(np)*dto2 
            py(np) = py(np) + vp(np)*dto2 
            pz(np) = pz(np) + wp(np)*dto2 
!
            uptilde(ijk) = up(np) 
            vptilde(ijk) = vp(np) 
            wptilde(ijk) = wp(np) 
         end do 
!
!
!     ******************************************************************
!
!     calculate new natural coordinates of mid-point
!
!     ******************************************************************
!
         call parlocat (ncellsp, ijkctmp, iphead, iwid, jwid, kwid, nsampl, &
            vrms, t_wall, ico, ijktmp2, ijktmp3, ijktmp4, rmaj, dz, divpix, &
            area_x, area_y, area_z, dto2, itdim, npart, ibar, jbar, kbar, mgeom&
            , cdlt, sdlt, wate, x, y, z, bxv, byv, bzv, xptilde, yptilde, &
            zptilde, uptilde, vptilde, wptilde, tsix, tsiy, tsiz, etax, etay, &
            etaz, nux, nuy, nuz, px, py, pz, up, vp, wp, pxi, peta, pzta, &
            nu_len, nu_comm) 
!
!
         call watev (ncellsp, ijkctmp, iphead, itdim, pxi, peta, pzta, wate) 
!
!
!     **************************************************
!
!
!     corrector step
!
!
!     *************************************************
!
         call trilinp (ncellsp, ijkctmp, itdim, iwid, jwid, kwid, iphead, pxi, &
            peta, pzta, wate, ex, ey, ez, expn, eypn, ezpn) 
!
         call trilinp (ncellsp, ijkctmp, itdim, iwid, jwid, kwid, iphead, pxi, &
            peta, pzta, wate, bxv, byv, bzv, bxpn, bypn, bzpn) 
!
         call trilinp (ncellsp, ijkctmp, itdim, iwid, jwid, kwid, iphead, pxi, &
            peta, pzta, wate, gradxb, gradyb, gradzb, uptilde, vptilde, wptilde&
            ) 
!
!
!
!     **************************************************************
!
!     solve the momentum equation
!
!     ************************************************************
!
!dir$ ivdep
!
         do n = 1, ncellsp 
!
            ijk = ijkctmp(n) 
            np = iphead(ijk) 
!
            is = ico(np) 
            qomdt = qom(is)*dto2 
!
            rqom = 1./qom(is) 
!
            omdtsq = qomdt**2*(bxpn(ijk)**2+bypn(ijk)**2+bzpn(ijk)**2) 
            denom = 1./(1. + omdtsq) 
!
!
!     solve the momentum equation
!
            ehatx = expn(ijk) - pmu(np)*uptilde(ijk)*rqom 
            ehaty = eypn(ijk) - pmu(np)*vptilde(ijk)*rqom 
            ehatz = ezpn(ijk) - pmu(np)*wptilde(ijk)*rqom 
!
            rbsq = 1./(bxpn(ijk)**2+bypn(ijk)**2+bzpn(ijk)**2) 
!
            udotb = (up(np)+qomdt*ehatx)*bxpn(ijk) + (vp(np)+qomdt*ehaty)*bypn(&
               ijk) + (wp(np)+qomdt*ehatz)*bzpn(ijk) 
!
            uptilde(ijk) = (udotb*bxpn(ijk)+(ehaty*bzpn(ijk)-ehatz*bypn(ijk)))*&
               rbsq 
!
            vptilde(ijk) = (udotb*bypn(ijk)+(ehatz*bxpn(ijk)-ehatx*bzpn(ijk)))*&
               rbsq 
!
            wptilde(ijk) = (udotb*bzpn(ijk)+(ehatx*bypn(ijk)-ehaty*bxpn(ijk)))*&
               rbsq 
!
!
         end do 
!
!dir$ ivdep
!
         do n = 1, ncellsp 
!
            ijk = ijkctmp(n) 
            np = iphead(ijk) 
!
!
            up(np) = 2.*uptilde(ijk) - up(np) 
            vp(np) = 2.*vptilde(ijk) - vp(np) 
            wp(np) = 2.*wptilde(ijk) - wp(np) 
!
            px(np) = px(np) + up(np)*dto2 
            py(np) = py(np) + vp(np)*dto2 
            pz(np) = pz(np) + wp(np)*dto2 
!
            uptilde(ijk) = up(np) 
            vptilde(ijk) = vp(np) 
            wptilde(ijk) = wp(np) 
!
         end do 
!
         call parlocat (ncellsp, ijkctmp, iphead, iwid, jwid, kwid, nsampl, &
            vrms, t_wall, ico, ijktmp2, ijktmp3, ijktmp4, rmaj, dz, divpix, &
            area_x, area_y, area_z, dt, itdim, npart, ibar, jbar, kbar, mgeom, &
            cdlt, sdlt, wate, x, y, z, bxv, byv, bzv, xptilde, yptilde, zptilde&
            , uptilde, vptilde, wptilde, tsix, tsiy, tsiz, etax, etay, etaz, &
            nux, nuy, nuz, px, py, pz, up, vp, wp, pxi, peta, pzta, nu_len, &
            nu_comm) 
!
!
         do n = 1, ncellsp 
            ijk = ijkctmp(n) 
            np = iphead(ijk) 
            if (np <= 0) cycle  
!
            inew = int(pxi(np)) 
            jnew = int(peta(np)) 
            knew = int(pzta(np)) 
!
            ijknew = (knew - 1)*kwid + (jnew - 1)*jwid + (inew - 1)*iwid + 1 
!
            iphead(ijk) = link(np) 
            link(np) = iphd2(ijknew) 
            iphd2(ijknew) = np 
!
         end do 
!
!     ********************************************************
!
!     a routine to locate a particle on a grid
!     ************************************************************************
!
!     all particles have been processed
!
!     ***********************************************************************
 
         if (.not.nomore) go to 1 
!
      endif 
!
!dir$ ivdep
!
      do n = 1, ncells 
         ijk = ijkcell(n) 
         iphead(ijk) = iphd2(ijk) 
         iphd2(ijk) = 0 
      end do 
!
 
!
      return  
!l    ------------------------------------------> return ----->>>
!
      end subroutine parmovgc 
